using Plots, LinearAlgebra
For each value of $z$, this polynomial has some number of roots that, of course, change as $z$ varies.
## Here α and β must have the same length
# For example: A-M 3rd order
# α = [1,-1,0]
# β = [5/12,8/12,-1/12]
# Leapfrog
# α = [1,-1,0]
# β = [0,2,0]
###################
ρ = (α,z) -> (z.^(length(α)-1:-1:0))'*α
σ = (β,z) -> (z.^(length(β)-1:-1:0))'*β
R = (α,β,z) -> ρ(α,z)/σ(β,z)
#5 (generic function with 1 method)
Recall that the "possible boundary" of the stability region $S$ is given by
$$ \frac{\rho(e^{i \theta})}{\sigma(e^{i\theta})}, \quad \theta \in [0,2 \pi). $$function find_roots(c) # supposing that the leading order coefficient is 1
# c contains the remaining coefficients
r = length(c)
A = zeros(Complex{Float64},r,r)
A[1,:] = -c
A[2:end,1:end-1] = A[2:end,1:end-1] + I # add identity matrix to lower-left block
return eigvals(A)
end
function check_condition(λ)
if maximum(abs.(λ)) > 1
return 0
else
for i = 1:length(λ)
if abs(λ[i]) ≈ 1. && sum(map(t -> λ[i] ≈ t,λ)) > 1
return 0
end
end
end
return 1
end
function compute_roots(α,β,z)
r = length(α)-1
c = α-z*β
if α[1]-z*β[1] ≈ 0.
λ = find_roots(c[3:end]/c[2])
else
λ = find_roots(c[2:end]/c[1]) # let's suppose that first and second coefficients don't vanish simultaneously
end
return λ
end
function root_condition(α,β,z)
return compute_roots(α,β,z) |> check_condition
end
function check_convergence(α,β) #supposing that α[1] = 1
r = length(α)-1
if sum(α) ≈ 0. && sum(β) ≈ sum(α.*(r:-1:0))
println("Method is consistent")
else
println("Method is inconsistent")
end
if root_condition(α,β,0.) == 1
println("Method is zero-stable")
else
println("Method is unstable")
end
end
function convergence_stability(α,β)
check_convergence(α,β)
θ = 0:0.01(1+rand()/10):2*π # random perturbation to avoid singularities
z = map(t -> R(α,β,exp(1im*t)),θ);
if abs(minimum(real(z)) - maximum(real(z))) < .1
xrange = [-4.,4.]
else
xrange = [minimum(real(z))-1,maximum(real(z))+1]
end
if abs(minimum(imag(z)) - maximum(imag(z))) < .1
yrange = [-4.,4.]
else
yrange = [minimum(imag(z))-1,maximum(imag(z))+1]
end
if xrange[1] < -10
xrange[1] = -4
end
if yrange[1] < -10
yrange[1] = -4
end
if xrange[2] > 10
xrange[2] = 4
end
if yrange[2] > 10
yrange[2] = 4
end
contourf(xrange[1]:0.01:xrange[2],yrange[1]:0.01(1+rand()/10):yrange[2],(x,y)-> root_condition(α,β,x+1im*y),colorbar=false)
plot!(real(z),imag(z),xlim=xrange,ylim=yrange,aspectratio=1,legend=false,lw=4,linecolor=:orange)
end
convergence_stability (generic function with 1 method)
Find the roots of
\begin{align*} \zeta^2 - \zeta = 0 \end{align*}find_roots([-1,0])
2-element Vector{ComplexF64}:
0.0 + 0.0im
1.0 + 0.0im
# A-M 3rd order
α = [1,-1,0]
β = [5/12,8/12,-1/12]
p = convergence_stability(α,β)
Method is consistent Method is zero-stable
savefig(p,"AM-3rd.pdf")
compute_roots(α,β,-2)
2-element Vector{ComplexF64}:
-0.40582741955797763 + 0.0im
0.2240092377397959 + 0.0im
# Leapfrog
α = [1,0,-1]
β = [0,2,0]
convergence_stability(α,β)
Method is consistent Method is zero-stable
# Forward Euler
α = [1,-1]
β = [0,1]
convergence_stability(α,β)
Method is consistent Method is zero-stable
# Backward Euler
α = [1,-1]
β = [1,0]
convergence_stability(α,β)
Method is consistent Method is zero-stable
# Trapezoid
α = [1,-1]
β = [.5,.5]
convergence_stability(α,β)
Method is consistent Method is zero-stable
# A-B four step
α = [1,-1,0,0,0]
β = [0,55,-59,37,-9]/24
convergence_stability(α,β)
Method is consistent Method is zero-stable
z = R(α,β,exp(1im*5))
0.40796674465342403 + 0.749343157823902im
abs.(compute_roots(α,β,z)) # roots of \pi(\xi;z)
4-element Vector{Float64}:
0.5341283948341674
1.000000000000001
0.37069069662575793
1.6159402431409078
z = R(α,β,exp(1im*4.252))
8.553795540246944e-5 + 0.43003577762135065im
compute_roots(α,β,z)
4-element Vector{ComplexF64}:
-0.4442966202159724 + 0.8958797426355074im
0.19049859439097486 - 0.4108913227581007im
0.3459692948015803 + 0.08832682135670405im
0.9080247555045485 + 0.4121834158148169im